

*** Replication .do File for "Diversity and Violence During Conflict Migration: The Troubles in Northern Ireland"
*** Adida, Brown, McCord and McLachlan

capture log close
log using "replication.txt", replace text
set more off

ssc install labsumm

use "data_final.dta", clear

********************* Pre-analysis variable construction

* Construct simple church diversity variable:
gen IV_ker0p5 = 1-(abs(cathker0p5-protker0p5))/(cathker0p5+protker0p5)
gen IV_ker1km = 1-(abs(cathker1-protker1))/(cathker1+protker1)
gen IV_ker1p5 = 1-(abs(cathker1p5-protker1p5))/(cathker1p5+protker1p5)
gen IV_ker2km = 1-(abs(cathker2-protker2))/(cathker2+protker2)
gen IV_ker3km = 1-(abs(cathker3-protker3))/(cathker3+protker3)
label var IV_ker0p5 "Church Diversity"
label var IV_ker1km "Church Diversity"
label var IV_ker1p5 "Church Diversity"
label var IV_ker2km "Church Diversity"
label var IV_ker3km "Church Diversity"

* Construct RQ version of church diversity
gen cathprop_ker1km = cathker1/(cathker1 + protker1)
gen cathprop_ker1p5 = cathker1p5/(cathker1p5 + protker1p5)
gen cathprop_ker2km = cathker2/(cathker2 + protker2)
gen cathprop_ker3km = cathker3/(cathker3 + protker3)
gen protprop_ker1km = 1 - cathprop_ker1km
gen protprop_ker1p5 = 1 - cathprop_ker1p5
gen protprop_ker2km = 1 - cathprop_ker2km
gen protprop_ker3km = 1 - cathprop_ker3km
gen rq_ker1km = 1 - (cathprop_ker1km * ( (0.5-cathprop_ker1km)/0.5)^2) - (protprop_ker1km * ( (0.5-protprop_ker1km)/0.5)^2)
gen rq_ker1p5 = 1 - (cathprop_ker1p5 * ( (0.5-cathprop_ker1p5)/0.5)^2) - (protprop_ker1p5 * ( (0.5-protprop_ker1p5)/0.5)^2)
gen rq_ker2km = 1 - (cathprop_ker2km * ( (0.5-cathprop_ker2km)/0.5)^2) - (protprop_ker2km * ( (0.5-protprop_ker2km)/0.5)^2)
gen rq_ker3km = 1 - (cathprop_ker3km * ( (0.5-cathprop_ker3km)/0.5)^2) - (protprop_ker3km * ( (0.5-protprop_ker3km)/0.5)^2)
label var rq_ker1km "Church Diversity (1 km bw)"
label var rq_ker1p5 "Church Diversity (1.5 km bw)"
label var rq_ker2km "Church Diversity (2 km bw)"
label var rq_ker3km "Church Diversity (3 km bw)"

* Herfindahl version of instrument:
gen hindex_ker1km = (cathker1/(cathker1 + protker1))^2 + (1-(cathker1/(cathker1 + protker1)))^2
gen hindex_ker1p5 = (cathker1p5/(cathker1p5 + protker1p5))^2 + (1-(cathker1p5/(cathker1p5 + protker1p5)))^2
gen hindex_ker2km = (cathker2/(cathker2 + protker2))^2 + (1-(cathker2/(cathker2 + protker2)))^2
gen hindex_ker3km = (cathker3/(cathker3 + protker3))^2 + (1-(cathker3/(cathker3 + protker3)))^2
label var hindex_ker1km "Church HHI"
label var hindex_ker1p5 "Church HHI"
label var hindex_ker2km "Church HHI"
label var hindex_ker3km "Church HHI"

* ELF version of instrument:
gen elf_ker1km =  1 - ((cathker1/(cathker1 + protker1))^2 + (1-(cathker1/(cathker1 + protker1)))^2)
gen elf_ker1p5 = 1 - ((cathker1p5/(cathker1p5 + protker1p5))^2 + (1-(cathker1p5/(cathker1p5 + protker1p5)))^2)
gen elf_ker2km = 1 - ((cathker2/(cathker2 + protker2))^2 + (1-(cathker2/(cathker2 + protker2)))^2)
gen elf_ker3km = 1 - ((cathker3/(cathker3 + protker3))^2 + (1-(cathker3/(cathker3 + protker3)))^2)
label var elf_ker1km "Church ELF"
label var elf_ker1p5 "Church ELF"
label var elf_ker2km "Church ELF"
label var elf_ker3km "Church ELF"

**** Prepare population variables
gen lpop = ln(population)
label var lpop "ln(Population)"
gen lpop2 = lpop^2
label var lpop2 "ln(Population) ^ 2"
gen lpop3 = lpop^3
gen lpop4 = lpop^4
gen nohotwater = v359_71/population
gen nobathshower = v360_71/population
gen noinsidetoilet = v361_71/population
gen notoilet = v362_71/population
label var nohotwater "Population without Hot Water"
label var nobathshower "Population without Bath or Shower in Home"
label var noinsidetoilet "Population without Inside Toilet"
label var notoilet "Population without Toilet"
label var victims "Victims in Grid Cell"

*** Create a PCA variable for the SES variables:
pca nohotwater nobathshower noinsidetoilet notoilet
predict pc1, score

*** Because of the decade-by-decade structure of data, there are now repeated observations per location
*** When running regressions with all victims, you'll get repeated observations, so use an index:
bysort boxcode: gen alldecades=1 if _n==1

gen isviolent = 1 if alldecades==1 & victims>0 & victims~=.
replace isviolent = 0 if isviolent==. & victims~=.

gen ldist2chca = ln( dist2chca )
label var ldist2chca "ln(Dist to nearest Catholic Church)"
gen ldist2chpr = ln( dist2chpr )
label var ldist2chpr "ln(Dist to nearest Protestant Church)"

/*** Likelihood ratio test to see if we need Neg Binomial
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop nohotwater nobathshower noinsidetoilet if alldecades==1
estimates store p1
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop nohotwater nobathshower noinsidetoilet if alldecades==1
estimates store nb1
lrtest p1 nb1, force
**/

******************** Table 2: Summary Statistics
labsumm victims rq_ker* dist2ch* pop nohotwater nobathshower noinsidetoilet if alldecades==1 

******* Figure 4: Polynomial plot of reduced form (victims on church diversity):
local BWs "1km 1p5 2km 3km"
foreach BWvar of local BWs {
	twoway (lpolyci victims rq_ker`BWvar', clcolor(maroon) ciplot(rline) clwidth(medthick) blpattern(dash) blcolor(gs6)) if alldecades==1, xtitle("Church Diversity") ytitle("Number of Victims") legend(off) 
	graph export "Fig4_`BWvar'.png", as(png) replace
}
**/

*** Table 3:   Negative Binomial estimates 
** Column 1: 1.5 km w/o SES
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop  if alldecades==1, vce(cluster ward93_id)
** Column 2: 2 km w/o SES
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 3: 3 km w/o SES
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 4: 1.5 km 
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 5: 2 km 
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 6: 3 km 
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/


*** Figure 5: Marginal Effects of Negative Binomial with SES
local BWs "1p5 2km 3km"
foreach BWvar of local BWs { 
	quietly nbreg victims c.rq_ker`BWvar'##c.rq_ker`BWvar' ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
	margins, dydx(rq_ker`BWvar') at(rq_ker`BWvar'=(0(0.1)1))
	marginsplot, title("") name(margins`BWvar', replace)
	marginsplot, title("") name(margins`BWvar', replace) addplot(hist rq_ker`BWvar', below width(0.1) yaxis(2) fcolor(gs12%30) lcolor(gs12%0) legend(off))
	graph export "Figure5_`BWvar'.png", as(png) replace
}
**/

*** Figure 6: Marginal plot on main regression by % Catholic density
graph drop _all
capture drop cathpropcat* protpropcat*
local BWs "1p5 2km 3km"
capture label define bins 10 "0-10" 20 "10-20" 30 "20-30" 40 "30-40" 50 "40-50" 60 "50-60" 70 "60-70" 80 "70-80" 90 "80-90" 100 "90-100"
foreach BWvar of local BWs {
	gen cathpropcat_`BWvar' = .
	replace cathpropcat_`BWvar' = 10 if (cathprop_ker`BWvar' <= 0.1)
	forvalues i = 0.2(0.1)1 {
		replace cathpropcat_`BWvar' = (`i'*100) if cathprop_ker`BWvar' <= `i' & cathpropcat_`BWvar'== .
	}
	capture replace cathpropcat_1p5 = 100 if cathprop_ker1p5 == 1
	label var cathpropcat_`BWvar' "Catholic Proportion of Church Density"
	label values cathpropcat_`BWvar' bins
	quietly nbreg victims c.rq_ker`BWvar'##c.rq_ker`BWvar' ldist2chca ldist2chpr if alldecades==1, vce(cluster ward93_id)
	margins, over(cathpropcat_`BWvar')
	marginsplot, title("") name(cathmargins_`BWvar', replace) xscale(range(10 100)) xlabel(10(10)100)
	graph export "Figure6_catholic_`BWvar'.png", as(png) replace
	
	gen protpropcat_`BWvar' = .
	replace protpropcat_`BWvar' = 10 if (protprop_ker`BWvar' <= 0.1)
	forvalues i = 0.2(0.1)1 {
		replace protpropcat_`BWvar' = (`i'*100) if protprop_ker`BWvar' <= `i' & protpropcat_`BWvar'== .
	}
	capture replace protpropcat_2km = 100 if protprop_ker2km == 1
	capture replace protpropcat_3km = 100 if protprop_ker3km == 1
	label values protpropcat_`BWvar' bins
	label var protpropcat_`BWvar' "Protestant Proportion of Church Density"
	quietly nbreg victims c.rq_ker`BWvar'##c.rq_ker`BWvar' ldist2chca ldist2chpr if alldecades==1, vce(cluster ward93_id)
	margins, over(protpropcat_`BWvar')
	marginsplot, title("") name(protmargins_`BWvar', replace) xscale(range(10 100)) xlabel(10(10)100)
	graph export "Figure6_protestant_`BWvar'.png", as(png) replace
	
}
**/


****** Figure A1: Histograms of church diversity variable by bandwidth:
twoway hist rq_ker1p5 if alldecades==1, name(hist_1p5km, replace) 
graph export "FigureA1_histogram1p5.png", as(png) replace
twoway hist rq_ker2km if alldecades==1, name(hist_2km, replace) 
graph export "FigureA1_histogram2.png", as(png) replace
twoway hist rq_ker3km if alldecades==1, name(hist_3km, replace) 
graph export "FigureA1_histogram3.png", as(png) replace
**/


*** Table A1: Main result using Poisson instead of Negative Binomial
** Column 1: 1.5 km w/o SES
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop  if alldecades==1, vce(cluster ward93_id)
** Column 2: 2 km w/o SES
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 3: 3 km w/o SES
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 4: 1.5 km 
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 5: 2 km 
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 6: 3 km 
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
poisson victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/


*** Figure A2: Marginal Effects of Poisson 
local BWs "1p5 2km 3km"
foreach BWvar of local BWs { // with SES
	quietly poisson victims c.rq_ker`BWvar'##c.rq_ker`BWvar' ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
	margins, dydx(rq_ker`BWvar') at(rq_ker`BWvar'=(0(0.1)1))
	marginsplot, title("") name(margins`BWvar', replace)
	marginsplot, title("") name(margins`BWvar', replace) addplot(hist rq_ker`BWvar', below width(0.1) yaxis(2) fcolor(gs12%30) lcolor(gs12%0) legend(off))
	graph export "FigureA2_margeffect_`BWvar'.png", as(png) replace
}
**/

*** Table A2: Results using raw census data
** Construct diversity variable from census
capture gen pctca = krig_pctca/100
capture gen pctpr = 1 - pctca
capture gen rq_diversity = 1 - (pctca * ( (0.5-pctca)/0.5)^2) - (pctpr * ( (0.5-pctpr)/0.5)^2)
capture gen rq_diversity_sq = rq_diversity^2
label var rq_diversity "Catholic/Protestant Diversity"
label var rq_diversity_sq "Catholic/Protestant Diversity Squared"
** Column 1: Endogenous regression without SES
nbreg victims rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 2: Endogenous regression with SES
nbreg victims rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/


** Table A3: Exclude Deaths later than 1980 or later than 1990 (Neg Binomial)
capture drop victims_pre1980 victims_pre1990
bysort boxcode: egen victims_pre1980 = sum(decade_victims) if decade==1960 | decade==1970
replace victims_pre1980=0 if victims_pre1980==.
bysort boxcode: egen victims_pre1990 = sum(decade_victims) if decade==1960 | decade==1970 | decade==1980
replace victims_pre1990=0 if victims_pre1990==.
** Column 1: 1.5 km All
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 2: 1.5 km pre-1990
nbreg victims_pre1990 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 3: 1.5 km pre-1980
nbreg victims_pre1980 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 4: 2 km All
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 5: 2 km pre-1990
nbreg victims_pre1990 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 6: 2 km pre-1980
nbreg victims_pre1980 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 7: 3 km All
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 8: 3 km pre-1990
nbreg victims_pre1990 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 9: 3 km pre-1980
nbreg victims_pre1980 churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
**/


*** Table A4: Two-panel analysis comparing church measure and census data measure, by decade.
** Church diversity variable
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker2km
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
** Construct diversity variable from census
capture gen pctca = krig_pctca/100
capture gen pctpr = 1 - pctca
capture gen rq_diversity = 1 - (pctca * ( (0.5-pctca)/0.5)^2) - (pctpr * ( (0.5-pctpr)/0.5)^2)
capture gen rq_diversity_sq = rq_diversity^2
label var rq_diversity "Catholic/Protestant Diversity"
label var rq_diversity_sq "Catholic/Protestant Diversity Squared"
** Construct victims by decade-by-decade
bysort boxcode: egen victims_1970s = sum(decade_victims) if decade==1970
replace victims_1970s=0 if victims_1970s==.
bysort boxcode: egen victims_1980s = sum(decade_victims) if decade==1980 
replace victims_1980s=0 if victims_1980s==.
bysort boxcode: egen victims_1990s = sum(decade_victims) if decade==1990
replace victims_1990s=0 if victims_1990s==.
bysort boxcode: egen victims_80s90s = sum(decade_victims) if decade==1980 | decade==1990
replace victims_80s90s=0 if victims_80s90s==.
** Column 1a: Endogenous regression for 1970s without SES
nbreg victims_1970s churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 2a: Endogenous regression for 1970s with SES
nbreg victims_1970s churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 3a: Endogenous regression for 1980s-1990s without SES
nbreg victims_80s90s churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 4a: Endogenous regression for 1980s-1990s with SES
nbreg victims_80s90s churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**
** Column 1b: Endogenous regression for 1970s without SES
nbreg victims_1970s rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 2b: Endogenous regression for 1970s with SES
nbreg victims_1970s rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 3b: Endogenous regression for 1980s-1990s without SES
nbreg victims_80s90s rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 4b: Endogenous regression for 1980s-1990s with SES
nbreg victims_80s90s rq_diversity rq_diversity_sq ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/


*** Table A5:  Main table excluding cells with "found" victims
gen foundcell = (victims > victims_exactloc)
** Column 1: 1.5 km w/o SES
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop  if alldecades==1 & foundcell==0, vce(cluster ward93_id)
** Column 2: 2 km w/o SES
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop  if alldecades==1 & foundcell==0, vce(cluster ward93_id)
** Column 3: 3 km w/o SES
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1 & foundcell==0, vce(cluster ward93_id)
** Column 4: 1.5 km w/SES
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1 & foundcell==0, vce(cluster ward93_id)
** Column 5: 2 km w/SES
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1 & foundcell==0, vce(cluster ward93_id)
** Column 6: 3 km w/SES
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1 & foundcell==0 , vce(cluster ward93_id)
**/


*** Table A6: Intra-communal vs Other Deaths
** Create dependent variable
capture gen victims_intracomm = victims_cath_on_cath + victims_prot_on_prot
capture gen victims_notintra = victims - victims_intracomm
** Columns 1-2: 1.5 km 
capture drop churchdiv churchdiv2
gen churchdiv = rq_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims_intracomm churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
nbreg victims_notintra churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 3-4: 2 km
replace churchdiv = rq_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims_intracomm churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
nbreg victims_notintra churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 5-6: 3 km
replace churchdiv = rq_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims_intracomm churchdiv churchdiv2 ldist2chca ldist2chpr lpop nohotwater nobathshower noinsidetoilet if alldecades==1, vce(cluster ward93_id)
nbreg victims_notintra churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/



******  Table A7: Neg Bin using ELF (with church diversity quadratic, dist2church in logs, controlling for SES):
** Column 1: 1.5 km noSES
capture drop churchdiv churchdiv2
gen churchdiv = elf_ker1p5
label var churchdiv "Church Diversity"
gen churchdiv2 = churchdiv^2
label var churchdiv2 "Church Diversity Squared"
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 2: 2 km noSES
replace churchdiv = elf_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 3: 3 km noSES
replace churchdiv = elf_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop if alldecades==1, vce(cluster ward93_id)
** Column 1: 1.5 km withSES
replace churchdiv = elf_ker1p5
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 2: 2 km withSES
replace churchdiv = elf_ker2km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
** Column 3: 3 km withSES
replace churchdiv = elf_ker3km
replace churchdiv2 = churchdiv^2
nbreg victims churchdiv churchdiv2 ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
**/


*** Figure A6: Marginal Effects of Negative Binomial using ELF
local BWs "1p5 2km 3km"
foreach BWvar of local BWs { // with SES
	quietly nbreg victims c.elf_ker`BWvar'##c.elf_ker`BWvar' ldist2chca ldist2chpr lpop pc1 if alldecades==1, vce(cluster ward93_id)
	margins, dydx(elf_ker`BWvar') at(elf_ker`BWvar'=(0(0.05)0.5))
	marginsplot, title("") name(margins`BWvar', replace)
	marginsplot, title("") name(margins`BWvar', replace) addplot(hist elf_ker`BWvar', below width(0.05) yaxis(2) fcolor(gs12%30) lcolor(gs12%0) legend(off))
	graph export "FigureA3_`BWvar'.png", as(png) replace
}
**/

